UnitVector<-function(n){matrix(rep(1,n),n,1)}

sympower <- function(x,pow) {
edecomp <- eigen(x)
roots <- edecomp$val
v <- edecomp$vec
d <- diag(roots^pow)
sympow <- v %*% d %*% t(v)
sympow
}

P <-function(x) {
y <- x %*% solve(t(x) %*% x) %*% t(x)
y
}

Q <- function(x) {
q <- diag(dim(x)[1]) - P(x)
q
}

MakeFullMatrix <- function(x) {L <- length(x)
p <- (sqrt(8*L + 1) -1)/2
y <- matrix(0,p,p)
count <- 0
for(i in 1:p ) {for (j in 1:i){count <- count+1
y[i,j] <- x[count]
y[j,i] <- x[count]
}
}
y
}

MakeExactData <- function(mu,sigma,n) {
p <- dim(sigma)[1]
x <- matrix(rnorm(n*p),n,p)
x <-x %*% sympower(cov(x),-1/2) %*% sympower(sigma,1/2)
diff <- matrix(0,1,p)
for(i in 1:p){
diff[1,i] <- mu[i]-mean(x[,i])
}
for(i in 1:n) {
  for(j in 1:p) {
    x[i,j] <- x[i,j] + diff[1,j]}
}
x
}

MultivariateNormalSample <-function(mu,sigma,n) {
p <- dim(sigma)[1]
x <- matrix(rnorm(n*p),n,p)
x <- x %*% sympower(sigma,1/2)
for(i in 1:n) {
  for(j in 1:p) {
    x[i,j] <- x[i,j] + mu[j]}
}
x
}

MakeFactorCorrelationMatrix <- function(F){
r <-F %*% t(F)
h <- diag(r)
u <- diag(1 - h)
r <- r+u
r
}

MakeFactorData <- function(F,n) {
r <-MakeFactorCorrelationMatrix(F)
mu<-matrix(0,1,dim(F)[1])
x <-MultivariateNormalSample(mu,r,n)
x
}

OrthogonalizeColumns <- function(x) {
x <- x %*% sympower(t(x) %*% x, -1/2)
x
}

OrthogonalRightUnit <- function(x) {
p <- dim(x)[1]
q <- dim(x)[2]
m <- q - p
decomp <-svd(x)
v <- decomp$u
d <- diag(decomp$d)
w1 <- decomp$v
w2 <- OrthogonalizeColumns(Q(w1) %*% matrix(rnorm(q*m),q,m))
S <- OrthogonalizeColumns( matrix(rnorm(m*m),m,m))
oru <- w1 %*% t(w1) + w2 %*% S %*% t(w2)
oru
}


MakeExactFactorData <- function(F, n){

r <-MakeFactorCorrelationMatrix(F)
mu<-matrix(0,1,dim(F)[1])
x <-MakeExactData(mu,r,n)
x
}